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ABSTRACT 


This study is concerned with development of a variational formulation and a 
procedure for computational solution for the shape optimal design of a two-dimensional 
linear elastic body, using a mixed finite element discretization. Shape optimal design is a 
problem that has interested many researchers in the last fifteen years. The subject has 
been surveyed in a number of review articles, see e.g. Haftka and Grandhi [11]. 

Zienkiewicz and Campbell [ see , e.g. 12 ] , were among the first to approach this 
problem using finite element methods. Subsequently this method has been applied widely 
to problems in shape optimal design[see, e.g.2,3,13- 15] , but only with mixed success. 
The finite element method based on the displacement formulation has two main 
disadvantages, (1) the increase of finite element error that results from mesh distortion 
during shape redesign, and (2) in some situations, a lack of sufficient precision in the 
prediction of stresses and strains at the boundary and internal nodes. There are some 
methods one can consider to overcome these difficulties. Some investigators have 
applied the Boundary Element Method[see, e.g. 16 - 18]. While the BEM has proved to be 
very useful and looks promising in certain applications of shape optimal design, for 
problems that require numerous evaluations of state variables in the domain (objective 


function = max xe QF(u,e), for example) the BEM loses some of its advantages, also at the 
current stage of development it lacks the generality provided by FEM in structural analysis. 


Within the FEM applications, the domain method [see, e.g. 19 ], where sensitivity 
expressions are defined in terms of domain integrals rather than boundary integrals 
(thereby avoiding the evaluation of state variables at the boundary ), provides for improved 
accuracy in the numerical calculation of sensitivities. Also recently, Haber proposed an 
Eulerian - Lagrangian formulation based on the mutual Reissner energy [see, e.g. 20], 
where the shape optimization problem can be formulated in an arbitrary initial domain as a 
means to overcome the difficulties inherited from shape redesign . 

In this work another approach is considered. With the development of automatic 
mesh generation and optimization techniques[see, e.g.7] the first of the cited disadvantages 
of FEM is avoided. Mixed finite element methods[see, e.g.8] that may provide for 
accurate computation of stresses and strains at the element nodes appears to be a natural 
approach to resolve the other difficulty. These considerations are brought together in the 
developments reported here, to demonstrate a more effective approach to the overall 
treatment of shape optimal design. 
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EXAMPLES OF SHAPE OPTIMIZATION PROBLEMS 


For simplicity let us consider plane linearly elastic structures to find the oprimal 
shape using the mixed finite element formulation together with an automatic mesh 
generation method based on the elliptic differential equations. Two model problems 
shown in the figure are solved by the present method to demonstrate its effectiveness. It is 
noted that the fillet optimal shape design problem (Model B) is a standard one, but is one of 
the most difficult problems because of sharp design change at the left edge of the design 
boundary. 
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DIFFICULTIES OF SHAPE OPTIMIZATION PROBLEMS 


The difficulty of the shape optimization problem is large design changes lead to 
significant changes of the corresponding finite element model of the structure during an 
optimization process. If the final, that is, optimal shape is known, it is possible to set the 
initial finite element model very close to the optimum. In this case, design change does not 
imply large change of finite element models, and then it is possible to avoid distorted finite 
elements which yield significant approximation errors and sometimes even negative values 
of the Jacobian of the isoparametric transformation. It is natural that the optimal shape is 
not known a priori, and then it is necessary to establish a shape optimization method to deal 
with large design change. Finite element approximation errors are strongly dependent 
upon the size and shape of finite elements. Errors are generally very large in regions 
where stresses are rapidly varying. In most shape optimization problems, stresses are 
varying rapidly at the end points of the design boundary where shape also changes rapidly. 
Furthermore, if shape change is large, this almost automatically yields distortion of finite 
elements, i.e. generation of unnecessary approximation errors. Figure (a) shows the 
distribution of finite element approximation error in a similar problem to Model A, while 
Figure (b) indicates a pathology in shape optimization. 



(b) Oscillation of the Design Boundary ( by. C. Fleury ) 
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FORMULATION OF THE DESIGN PROBLEM 

The optimal shape design problem is defined by minimizing the maximum value of 

a local criterion F(u,e) stated by the displacement u and the linearized strain tensor e, 
subject to the state equations which represent equilibrium, constitutive relation, and 
boundary conditions. Here D represents the design variable to describe the shape of the 
design boundary. 


Mino MaxxeQ F (u(x),e(x)) 
subject to 

the resource constraint J dQ - A < 0 

Q 

the equilibrium equations div Z + f = 0 , E = E T in Q 
the strain - displacement relations e = j(Vu + uV) in Q. 

the stress - strain relation X = E : e in Q 

die traction boundary condition nZ = t on T t 

the displacement boundary condition u = 0 on T u 
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MIXED FINITE ELEMENT FORMULATION 

If the displacement method is applied in finite element analysis, strain and stress 
components are computed at each Gaussian integration point to form element stiffness 
matrices and load vectors. However, values of stress components must be obtained at the 
nodes on the design boundary for the shape optimization problems. Thus an extrapolation 
method must be introduced to obtain nodal values of stress components. For example, if 
the least squares method is applied to obtain nodal values, it cannot provide sufficiently 
accurate values in the region that stress gradient is high. Furthermore, it becomes very 
inaccurate if distorted irregular finite elements exist in a finite element model. To avoid 
such problems, we here apply a mixed formulation that computes nodal values of 
displacements, strains, and stresses, directly without applying an extrapolation method. 
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QUAD4 

8 degrees of freedom per node { u x ,u y ,e X x,eyy,Txy,cJxx,^yy ,a xy } T 

4 

u x = £ a=1 UxaNa(^dl)> etc 

4 

£xx = £ a= i ExxaNa^/nX etc 

4 

Gxx = ^a=l a xxaNa(£/n)> etc 
Na(^ri) = \ (l+^aW+tlcdl), a=l,...,4 
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OPTIMALITY CONDITIONS 


Transferring the original shape optimization problem to the upper-bound 
formulation MinD P subject to the above constraints and the additional one F(u,e) - P < 0 
in Q, the Lagrange multiplier method implies the necessary condition for optimality as 
shown in below. 


Equilibrium Equations 


f(5<V8v+6vV):Z-5v.f}dO - }5v.tdr 
3 r ‘ 

+ J[Se:(Ee-E) + 8T:{|(Vu+uV)-e}]dQ = 0, u=0 on r u . 
Q 


Adjoint Problems 


J{l<Vv+vV):5E + X(||: 8e + g|.8u)}dft 


n 


+ 


Je:( E: 


8e - 8E) + T:(i<V8u+8uV) - 8e)} dQ = 0 
2 


n 


Optimality Condition due to Variation by Shape Change 

f (Z:e+ A)(0.n) dT = 0, A(fdQ-A)=0, 0<A, fdQ-A<0 

r d n q 

Normality Condition 

J fcdQ = 1, X(F(u,e) - p) = 0, 0<X, F(u,e) - P<0, in Q 
Q 
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SOLUTION PROCEDURE : OPTIMIZATION ALGORITHM 
Solution procedure for the shape optimization is described in the following flow 

chart. 
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MESH GENERATION METHOD : ELLIPTIC MESH GENERATOR 


At each design iteration, a new finite element mesh is reconstructed using the elliptic 
mesh generation method. In most of shape optimization practice, finite element meshes are 
not regenerated, but are modified during the optimization process. Because of the method 
applied for modification of the initial finite element mesh, unnecessary element distortion 
is, in general, generated in each iteration. Element distortion can become so large it can 
actually destroy accuracy of the finite element approximations, and then yield 
unsatisfactory results in shape optimization. To avoid this difficulty, it is better to 
regenerate a finite element mesh at each design stage despite of expense required. Here the 
elliptic mesh generation method is used that generates almost orthogonal meshes using only 
the data of the boundary of the domain occupied by a structure. If the boundary is 
represented by a set of spline curves, it is possible to represent the design boundary by 

several spline functions defined by the location 0 a , a=l,...,a max , of the so-called control 
points without loss of generality. 


Elliptic Mesh Generation Method 

V(DV*) — PfotO and V(DVq) = -Q£,ri) 

3 3 

where V = Uj— + (x,y) are the physical coordinates, and (^,r|) are the 

mesh coordinates. 


I I I I I t I 

• - 1 - f -I 1- 

I I I I I I I 

— i “ "if - r "i — i — k “ r “ 

i i i i i i i 

"vrrrvrr 

- j - i. .i — i. j - i . 

i i i i i i i 

- H - i - f -I- -I- 

I I I I I I I 

" t - r "r "i“ t r “ 

i i i i i t i 

“ ‘« " T T T "■ ’ V “ 

O - 1 — ‘-’-O— ^ 


Mesh Coordinates 



O Control Points for Mesh 
on the Design Boundary 
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EXAMPLE: MODEL A 


Model Problem A is to find the optimally shaped hole in a biaxially loaded linearly 
elastic thin plate. Starting from a rectangular hole, the optimal shape of the hole is 
obtained for the loading condition which yields an elliptic hole as the optimum. Applying 
the symmetry condition, only a quarter part of the plate is discretized by 1 10 QUAD4 
elements together with 1 1 control points on the design boundary for shape optimization. 
The minimum admissible area of the hole prescribed by the user is restricted to A=0.16m 2 . 
Iteration history and the optimal shape are given in the figure. The first three redesigns 
rapidly reduce the maximum value of the von Mises equivalent stress from 6.6 to 3.3, 
while the optimum obtained after 23 iterations is about 2.8. 


Objective Function 



Iteration 23 

Iteration 1 
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EXAMPLE: MODEL B 


For finding the optimal shape for a linearly elastic fillet, the domain is discretized by 
126 QUAD4 elements together with 14 control points on the design boundary. In this 
case, the objective function F is the von Mises equivalent stress. Assuming that the 
maximum allowable area of the fillet is restricted to A=1.135 x 10' 2 m 2 , an optimal shape is 
obtained. The maximum value of the equivalent stress becomes very stable after 15 design 
iterations, and its minimum is achieved at the 24th iteration. This means that optimization 
and finite element remeshing is performed 24 times. First 10 iterations give rapid 
reduction of the maximum value of the stress. The maximum von Mises stress in the final 
fillet obtained is less than half of the maximum stress of the initial structure. 
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Iteration 11 


Iteration 24 


FINAL REMARKS 


The development presented introduces a general form of mixed formulation for the 
optimal shape design problem, and the associated optimality conditions are easily obtained 
without resorting to highly elaborate mathematical developments. Also the physical 
significance of the adjoint problem comes out to be clearly defined with this formulation. 

In the examples presented, an elliptical automatic mesh generator assuring an 
orthogonal finite element mesh at the domain boundary [ see, e.g. 7] was used at each 
shape redesign . Although this procedure might seem to be computationally a very 
expensive procedure, actually it guarantees a good accuracy for the discrete model with an 
increase on computational time of less than 5 % of the actual time required for the finite 
element analysis. 

The numerical examples presented demonstrates the stability of the procedure. 
Problems commonly encountered in shape optimization arising from the development of 
instabilities in the design boundary definition were largely avoided. As is to be expected, 
however, this improvement is accomplished at the expense of the increase in cost of 
computation as compared to the simple displacement formulation. 


SHAPE OP TIMIZATION 

1. Mixed Finite Element Methods for Analysis 

2. Automatic Remeshing Scheme by 
Elliptic Mesh Generation Methods 

3. Optimality Conditions are Obtained by the Upper Bound 
Method 

4. Demonstration by Numerical Examples 
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